#include <iostream>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <limits.h>
#include <float.h>
#include <math.h>
#include <cmath>

#define epochs_num 1000
#define learning_rate_att 0.01
#define learning_rate_rep 0.1
#define L 5

using namespace std;

double x[1000][1000];
int similarity[1000][1000];
int label[1000];
double weight[10][10][10];
double y_p[10][10];
double y_q[10][10];
double y_r[10][10];
double y_s[10][10];
double x_p[10][10];
double x_q[10][10];
double x_r[10][10];
double x_s[10][10];
int stop_layer = 1;
int att_flag = 1;
int rep_flag = 1;

FILE* train_file;

int read_input(FILE* f, double x[1000][1000]);
void read_similarity(FILE *f, int data_n, int similarity[1000][1000]);
void print_input(int data_n, double x[1000][1000]);
void print_similarity(int data_n, int similarity[1000][1000]);
void init_weight(int layer, int neuron_n[], double weight[10][10][10]);
void print_weight(int layer, int neuron_n[], double weight[10][10][10]);
void calc_layer_m_output(int m, int neuron_n[], double input[], double weight[10][10][10], double output[]);
double calc_distance(int n, double x[], double y[]);
void assign_y_p_y_q(int n, double aa1[], double aa2[], double y1[], double y2[]);
void update_weight(int m, int neuron_n[], double w[10][10]);

int main()
{
    int neuron_n[7] = {2,5,5,5,5,5,1};//input layer: 2 neurons output layer: 1 neuron
    int data_n = 0;
    char* pt_filename = "hw2pt.dat";
    char* class_filename = "hw2class.dat";
    char* output_filename = "hw2output.dat";
    FILE* hw2pt = fopen(pt_filename,"r");
    FILE* hw2class = fopen(class_filename,"r");
    //FILE* hw2output = fopen(output_filename,"w");
    train_file = fopen("train_file.dat","w");

    /**read data**/
    data_n = read_input(hw2pt, x);
    read_similarity(hw2class,data_n,similarity);

    /**init weight**/
    init_weight(L,neuron_n,weight);

    int m;
    int epochs;
    int i,j;

    for(m=1;m<=L;m++)
    {
        att_flag = 1;
        rep_flag = 1;
        double max_dist = 0;
        double min_dist = 1000;
        for(epochs=0;epochs<epochs_num;epochs++)
        {
            for(i = 0; i < 100; i++)
            {
                double y1[10];
                calc_layer_m_output(m,neuron_n,x[i],weight,y1);
                for(j = i + 1; j < 100; j++)
                {
                    double y2[10];
                    calc_layer_m_output(m,neuron_n,x[j],weight,y2);
                    double dist = 0;
                    dist = calc_distance(neuron_n[m], y1, y2);

                    if(similarity[i][j]==1&&i!=j)
                    {
                        if(dist>max_dist&&i!=j)
                        {
                            max_dist = dist;
                            assign_y_p_y_q(neuron_n[m],y_p[m],y_q[m],y1,y2);
                            memset(y1,0,sizeof(y1));
                            memset(y2,0,sizeof(y2));
                             //x_p x_q
                            int bottom_n = m-1;
                            calc_layer_m_output(bottom_n,neuron_n,x[i],weight,y1);
                            calc_layer_m_output(bottom_n,neuron_n,x[j],weight,y2);
                            assign_y_p_y_q(neuron_n[bottom_n],x_p[bottom_n],x_q[bottom_n],y1,y2);
                            //printf("max: i=%d j=%d\n",i,j);
                        }
                    }
                    else if(similarity[i][j]==0&&i!=j)
                    {
                        int n;
                        if(dist<min_dist)
                        {
                            min_dist = dist;
                            assign_y_p_y_q(neuron_n[m],y_r[m],y_s[m],y1,y2);
                            memset(y1,0,sizeof(y1));
                            memset(y2,0,sizeof(y2));
                            int bottom_n = m-1;
                            calc_layer_m_output(bottom_n,neuron_n,x[i],weight,y1);
                            calc_layer_m_output(bottom_n,neuron_n,x[j],weight,y2);
                            assign_y_p_y_q(neuron_n[bottom_n],x_r[bottom_n],x_s[bottom_n],y1,y2);
                            //printf("min: i=%d j=%d\n",i,j);
                        }
                    }
                }
            }
            if(min_dist >= (4*neuron_n[m]-1)){rep_flag = 0;}
            else if(min_dist <= 1){rep_flag=1;}
            if(max_dist <= 1){att_flag = 0;}
            else if(max_dist >= 3*neuron_n[m]){att_flag=1;}
            if(att_flag==0&&rep_flag==0){break;}
            update_weight(m, neuron_n, weight[m]);
            //printf("min_dist = %.2lf max_dist = %.2lf\n",min_dist,max_dist);
        }
        if(min_dist >= (4*neuron_n[m]-1)&&max_dist<=1){stop_layer=m;break;}
    }

    double ans[1000][10];
    int class_num = 0;
    for(i=0;i<100;i++){label[i]=-1;}
    for(i=0;i<100;i++)
    {
        for(j=i+1;j<100;j++)
        {
            if(similarity[i][j]==1&&label[i]==-1)
            {
                label[i]  = class_num;
                label[j]  = class_num;
                class_num++;
            }
            else if(similarity[i][j]==1&&label[i]!=-1)
            {
                label[j]  = label[i];
            }
        }
    }

    /*for(i=0;i<100;i++)
    {
        fprintf(hw2output,"%lf %lf %d\n",x[i][0],x[i][1],label[i]);
    }*/
    //FILE* matlab_output = fopen("matlab_output.dat","w");
    for(i=0;i<100;i++)
    {
        /*
        fprintf(train_file,"%lf %lf\n",x[i][0],x[i][1],label[i]);
        for(j=1;j<=L;j++)
        {
            calc_layer_m_output(j,neuron_n,x[i],weight,ans[i]);
            fprintf(train_file,"%lf %lf %lf %lf %lf %d\n",ans[i][0],ans[i][1],ans[i][2],ans[i][3],ans[i][4],label[i]);
        }
        fprintf(train_file,"\n");
        */
        calc_layer_m_output(stop_layer,neuron_n,x[i],weight,ans[i]);
        fprintf(train_file,"%lf %lf %lf %lf %lf %d\n",ans[i][0],ans[i][1],ans[i][2],ans[i][3],ans[i][4],label[i]);
    }
    print_weight(L,neuron_n,weight);

    return 0;
}

int read_input(FILE* f, double x[1000][1000])
{
    int n=0;
    while(fscanf(f,"%lf %lf",&x[n][0],&x[n][1])!=EOF)
    {
        n++;
    }
    return n;
}

void print_input(int data_n, double x[1000][1000])
{
    int i = 0;
    for(i=0;i<data_n;i++)
    {
        printf("%lf %lf\n",x[i][0],x[i][1]);
    }
}

void read_similarity(FILE *f, int data_n, int similarity[1000][1000])
{
    int i = 0;
    int j = 0;
    for(i=0;i<data_n;i++)
    {
        for(j=0;j<data_n;j++)
        {
            fscanf(f,"%d",&similarity[i][j]);
        }
    }
}

void print_similarity(int data_n, int similarity[1000][1000])
{
    int i = 0;
    int j = 0;
    for(i=0;i<data_n;i++)
    {
        for(j=0;j<data_n;j++)
        {
           printf("%d ",similarity[i][j]);
        }
        printf("\n");
    }
}

void init_weight(int layer, int neuron_n[], double weight[10][10][10])
{
    srand(time(NULL));
    int i;
    int j;
    int k;
    for(i=1;i<=L;i++)
    {
        for(j=0;j<neuron_n[i];j++)
        {
            for(k=0;k<neuron_n[i-1]+1;k++)
            {
                weight[i][j][k] = (2*((rand()%2)+1)-1);
            }
        }
    }
}

void print_weight(int layer, int neuron_n[], double weight[10][10][10])
{
    int i;
    int j;
    int k;
    for(i=1;i<=L;i++)
    {
        printf("layer %d:\n",i);
        for(j=0;j<neuron_n[i];j++)
        {
            printf("neuron %d: ",j);
            for(k=0;k<neuron_n[i-1]+1;k++)
            {
                printf("%.2lf ",weight[i][j][k]);
            }
            printf("\n");
        }
        printf("\n");
    }
}

void calc_layer_m_output(int m, int neuron_n[], double input[], double weight[10][10][10], double output[])
{
    int i,j,k;
    int n;

    double prev_layer_m_output[10];
    double curr_layer_m_output[10];

    prev_layer_m_output[0] = input[0];
    prev_layer_m_output[1] = input[1];
    prev_layer_m_output[2] = 1;

    if(m==0)
    {
        output[0] = input[0];
        output[1] = input[1];
        output[2] = 1;
        return;
    }

    for(i = 1; i <= m; i++)
    {
        for(j = 0; j < neuron_n[i]; j++)
        {
            double total = 0;
            for(k = 0; k < neuron_n[i-1] + 1; k++)
            {
                total = total + prev_layer_m_output[k]*weight[i][j][k];
            }
            curr_layer_m_output[j] = (2/(1+exp(-total)))-1;
        }
        for(n=0;n<neuron_n[i];n++){prev_layer_m_output[n] = curr_layer_m_output[n];}
        prev_layer_m_output[neuron_n[i]] = 1;
    }

    //output for the m-th layer's neuron that always output 1
    for(i=0;i<neuron_n[m];i++){output[i] = curr_layer_m_output[i];}
    output[neuron_n[m]] = 1;
}

double calc_distance(int n, double x[], double y[])
{
    int i,j;
    double dist = 0;
    for(i = 0; i < n; i++)
    {
        dist = dist + (x[i] - y[i])*(x[i] - y[i]);
    }
    return dist;
}

void assign_y_p_y_q(int n, double aa1[], double aa2[], double y1[], double y2[])
{
    int i;
    for(i = 0; i < n; i++)
    {
        aa1[i] = y1[i];
        aa2[i] = y2[i];
    }
}

void update_weight(int m, int neuron_n[], double w[10][10])
{
    int i,j;
    double der_E_att[10][10];
    double der_E_rep[10][10];

    for(i = 0; i < neuron_n[m]; i++)
    {
        for(j = 0; j < neuron_n[m-1] + 1; j++)
        {
            double const1 = 1 - (y_p[m][i]*y_p[m][i]);
            double const2 = 1 - (y_q[m][i]*y_q[m][i]);
            double const3 = 1 - (y_r[m][i]*y_r[m][i]);
            double const4 = 1 - (y_s[m][i]*y_s[m][i]);
            //printf("%.2lf %.2lf %.2lf %.2lf\n",const1,const2,const3,const4);
            if(j==neuron_n[m-1])
            {
                double yp_minus_yq = y_p[m][i] - y_q[m][i];
                double yr_minus_ys = y_r[m][i] - y_s[m][i];

                //derivative of E_att
                der_E_att[i][j] = yp_minus_yq*const1*(-1) - yp_minus_yq*const2*(-1);

                //derivative of E_rep
                der_E_rep[i][j] = -yr_minus_ys*const3*(-1)+yr_minus_ys*const4*(-1);
            }
            else
            {
                double yp_minus_yq = y_p[m][i] - y_q[m][i];
                double yr_minus_ys = y_r[m][i] - y_s[m][i];
                //derivative of E_att
                der_E_att[i][j] = yp_minus_yq*const1*x_p[m-1][j] - yp_minus_yq*const2*x_q[m-1][j];

                //derivative of E_rep
                der_E_rep[i][j] = -yr_minus_ys*const3*x_r[m-1][j] + yr_minus_ys*const4*x_s[m-1][j];
            }
        }
        //printf("%.2lf %.2lf\n",y_r[m][i],y_s[m][i]);
    }

    for(i = 0; i < neuron_n[m]; i++)
    {
        for(j = 0; j < neuron_n[m-1] + 1; j++)
        {
            double adjust = 0;
            if(att_flag==1)
            {
                adjust = adjust + learning_rate_att*der_E_att[i][j];
            }
            if(rep_flag==1)
            {
                adjust = adjust + learning_rate_rep*der_E_rep[i][j];
            }
            w[i][j] = w[i][j] - adjust;
        }
    }
}
